function [ests,params,gamma,kappa,sigma,b]=make_ests(theta,param_sims,error_sims,b_Omega,w_Omega)

b_mu=exp(theta(1))/(1+exp(theta(1)));
b_alpha=theta(2);
b_beta=theta(3);
b_sigma=exp(theta(4));
b_gamma_mn=theta(5);
b_gamma_sd=exp(theta(6));
w_mu=exp(theta(7))/(1+exp(theta(7)));
w_alpha=theta(8);
w_beta=theta(9);
w_sigma=exp(theta(10));
w_gamma_mn=theta(11);
w_gamma_sd=exp(theta(12));
corr_gamma=(exp(theta(13))-1)/(exp(theta(13))+1);
corr_sigma=(exp(theta(14))-1)/(exp(theta(14))+1);
L=length(param_sims)/2;

gamma_V=[b_gamma_sd^2,corr_gamma*b_gamma_sd*w_gamma_sd;corr_gamma*b_gamma_sd*w_gamma_sd,w_gamma_sd^2];
sigma_V=[b_sigma^2,corr_sigma*b_sigma*w_sigma;corr_sigma*b_sigma*w_sigma,w_sigma^2];

gamma=[b_gamma_mn,w_gamma_mn]+param_sims(1:L,:)*chol(gamma_V); 
kappa=[b_alpha,w_alpha]+repmat([b_beta,w_beta],L,1).*gamma+param_sims(L+1:end,:)*chol(sigma_V);

sigma=exp(kappa);
pi=zeros(L,2);
b=zeros(L,2);
for l=1:L
    pi(l,1)=signalcdf(gamma(l,1),b_mu,sigma(l,1));
    pi(l,2)=signalcdf(gamma(l,2),w_mu,sigma(l,2));
    b(l,1)=1/(1+((1-b_mu)/b_mu)*exp((0.5-gamma(l,1))/sigma(l,1)^2));
    b(l,2)=1/(1+((1-w_mu)/w_mu)*exp((0.5-gamma(l,2))/sigma(l,2)^2));
end

type1=zeros(L,2);
type2=zeros(L,2);
delta=zeros(L,2);
for l=1:L
    type1(l,1)=normcdf(-0.5/sigma(l,1)-sigma(l,1)*log(((1-b(l,1))*b_mu)/(b(l,1)*(1-b_mu))));
    type1(l,2)=normcdf(-0.5/sigma(l,2)-sigma(l,2)*log(((1-b(l,2))*w_mu)/(b(l,2)*(1-w_mu))));
    type2(l,1)=normcdf(-0.5/sigma(l,1)+sigma(l,1)*log(((1-b(l,1))*b_mu)/(b(l,1)*(1-b_mu))));
    type2(l,2)=normcdf(-0.5/sigma(l,2)+sigma(l,2)*log(((1-b(l,2))*w_mu)/(b(l,2)*(1-w_mu))));
    delta(l,1)=type1(l,1)*b_mu/pi(l,1);
    delta(l,2)=type1(l,2)*w_mu/pi(l,2);
end

params=[pi;delta];
ests=[params(:,1)+chol(b_Omega)'*error_sims(:,1),params(:,2)+chol(w_Omega)'*error_sims(:,2)];

